# import data
data <- read.csv("{TEMP_DIR}/data.csv")

# set row names
rownames(data) <- data$Sample

# delete column which held the row names
data$Sample <- NULL

# try to do factor analysis without obtaining covmat
try(fa <- factanal(x=data, factors={NUMBER_OF_FACTORS}, scores = "{SCORES}", rotation = "{ROTATION}"), TRUE)

if (!exists("fa")) {

    # load factor.model.stat function
    # "factor.model.stat" is an R (and S-PLUS) function that computes a variance matrix based on a factor model. This can estimate a variance model even when there are more variables than observations.
    # Source: http://www.burns-stat.com/
    "factor.model.stat" <-
    function (x, weight=seq(1, 3, length=nobs), output="full", center=TRUE, 
              frac.var=.5, iter.max=1, nfac.miss=1, full.min=20, reg.min=40, 
              sd.min=20, quan.sd=.90, tol=1e-3, zero.load=FALSE, 
              range.factors=c(0, Inf))
    {
        fun.copyright <- "Placed in the public domain in 2006 by Burns Statistics"
        fun.version <- "factor.model.stat 006"
        
        subfun.ssd <- function(z, weight, sd.min) {
            nas <- is.na(z)
            if(any(nas)) {
                if(sum(!nas) < sd.min) return(NA)
                    sum(weight[!nas] * z[!nas]^2) / sum(weight[!nas])
            } else {
                sum(weight * z^2)
            }
        }
            
    #
    # start of main function
    #
            
        if(is.data.frame(x)) {
            x <- as.matrix(x)
        } else {
            if(!is.matrix(x)) stop("x needs to be a matrix")
        }
            if(!is.numeric(x)) stop("x needs to be numeric")
                x[!is.finite(x)] <- NA
                
                xna <- is.na(x)
                allmis <- rowSums(xna) == ncol(x)
                if(any(allmis)) {
                    x <- x[!allmis, , drop=FALSE]
                    xna <- is.na(x)
                }
                num.mis <- colSums(xna)
                
                if(any(num.mis > 0)) {
                    if(sum(num.mis == 0) < full.min) 
                        stop("not enough columns without missing values")
                            if(!length(dimnames(x)[[2]])) 
                                stop("x needs column names when missing values exist")
                                    
                                    max.miss <- max(num.mis)
                                    lnfm <- length(nfac.miss)
                                    if(lnfm == 0) stop("nfac.miss must have positive length")
                                        nfac.miss <- round(nfac.miss)
                                        if(any(nfac.miss < 0)) 
                                            stop("negative values in nfac.miss")
                                                if(lnfm < max.miss) {
                                                    nfac.miss <- c(nfac.miss, rep(nfac.miss[lnfm],
                                                                                  max.miss - lnfm))
                                                }
                }
                                    
                                    if(length(output) != 1) stop("output needs to be a single string")
                                        output.menu <- c("full", "factor")
                                        output.num <- pmatch(output, output.menu, nomatch=0)
                                        if(output.num == 0) stop("unknown or ambiguous value for output")
                                            output <- output.menu[output.num]
                                            
                                            nassets <- ncol(x)
                                            nobs <- nrow(x)
                                            if(length(weight) != nobs) {
                                                if(length(weight) == nobs + sum(allmis)) {
                                                    weight <- weight[!allmis]
                                                } else {
                                                    stop("weight vector is the wrong length")
                                                }
                                            }
                                            if(any(weight < 0)) stop("negative weights not allowed")
                                                weight <- weight / sum(weight)
                                                
                                                if(is.logical(center)) {
                                                    if(center) {
                                                        center <- colMeans(x, na.rm=TRUE)
                                                    } else {
                                                        center <- rep(0, nassets)
                                                    }
                                                } else if(length(center) != nassets) stop("center is the wrong length")
                                                    x <- sweep(x, 2, center, "-")
                                                    
                                                    sdev <- sqrt(apply(x, 2, subfun.ssd, weight=weight, sd.min=sd.min))
                                                    if(any(sdev <= 0, na.rm=TRUE)) {
                                                        stop(paste(sum(sdev <= 0, na.rm=TRUE), 
                                                                   "asset(s) with constant returns"))
                                                    }
                                                    if(any(is.na(sdev))) {
                                                        sdev[is.na(sdev)] <- quantile(sdev, quan.sd, na.rm=TRUE)
                                                    }
                                                    x <- scale(x, scale=sdev, center=FALSE)
                                                    
                                                    xw <- sqrt(weight) * x
                                                    fullxw <- xw[, num.mis == 0, drop=FALSE]
                                                    fx.svd <- svd(fullxw, nu=0)
                                                    cumvar <- cumsum(fx.svd$d^2) / sum(fx.svd$d^2)
                                                    nfac <- sum(cumvar < frac.var) + 1
                                                    if(nfac > max(range.factors)) nfac <- max(range.factors)
                                                        if(nfac < min(range.factors)) nfac <- min(range.factors)
                                                            if(nfac > length(cumvar)) nfac <- length(cumvar)
                                                                fseq <- 1:nfac
                                                                loadings <- scale(fx.svd$v[, fseq, drop=FALSE], scale=fx.svd$d[fseq], 
                                                                                  center=FALSE)
                                                                
                                                                if(iter.max > 0) {
                                                                    cormat <- t(fullxw) %*% fullxw
                                                                    uniqueness <- 1 - rowSums(loadings^2)
                                                                    uniqueness[uniqueness < 0] <- 0
                                                                    uniqueness[uniqueness > 1] <- 1
                                                                    start <- uniqueness
                                                                    converged <- FALSE
                                                                    for(i in 1:iter.max) {
                                                                        cor.red <- cormat
                                                                        diag(cor.red) <- diag(cor.red) - uniqueness
                                                                        t.eig <- eigen(cor.red)
                                                                        t.val <- t.eig$value[fseq]
                                                                        t.val[t.val < 0] <- 0
                                                                        loadings <- scale(t.eig$vector[, fseq, drop=FALSE], 
                                                                                          center=FALSE, scale=1/sqrt(t.val))
                                                                        uniqueness <- 1 - rowSums(loadings^2)
                                                                        uniqueness[uniqueness < 0] <- 0
                                                                        uniqueness[uniqueness > 1] <- 1
                                                                        if(all(abs(uniqueness - start) < tol)) {
                                                                            converged <- TRUE
                                                                            break
                                                                        }
                                                                        start <- uniqueness
                                                                    }
                                                                }
                                                                dimnames(loadings) <- list(dimnames(fullxw)[[2]], NULL)
                                                                
                                                                if(any(num.mis)) {
    # calculate loadings for columns with NAs
                                                                    floadings <- loadings
                                                                    if(zero.load) {
                                                                        loadings <- array(0, c(nassets, nfac))
                                                                    } else {
                                                                        meanload <- colMeans(floadings)
                                                                        loadings <- t(array(meanload, c(nfac, nassets)))
                                                                    }
                                                                    dimnames(loadings) <- list(dimnames(x)[[2]], NULL)
                                                                    loadings[dimnames(floadings)[[1]], ] <- floadings
                                                                    scores <- fullxw %*% floadings
                                                                    dsquare <- fx.svd$d[1:nfac]^2
                                                                    nfac.miss[nfac.miss > nfac] <- nfac
                                                                    
                                                                    for(i in (1:nassets)[num.mis > 0 & nobs - num.mis > reg.min]) {
                                                                        t.nfac <- nfac.miss[ num.mis[i] ]
                                                                        if(t.nfac == 0) next
                                                                            t.okay <- !is.na(xw[, i])
                                                                            t.seq <- 1:t.nfac
                                                                            t.load <- lsfit(xw[t.okay, i], scores[t.okay, t.seq], 
                                                                                            intercept=FALSE)$coef / dsquare[t.seq]
                                                                            loadings[i, t.seq] <- t.load
                                                                            NULL
                                                                    }
                                                                }
                                                                    
                                                                    comm <- rowSums(loadings^2)
                                                                    if(any(comm > 1)) {
    # adjust loadings where communalities too large
                                                                        toobig <- comm > 1
                                                                        loadings[toobig,] <- loadings[toobig,] / sqrt(comm[toobig])
                                                                        comm[toobig] <- 1
                                                                    }
                                                                    
                                                                    switch(output, 
                                                                           full= {
                                                                               ans <- loadings %*% t(loadings)
                                                                               diag(ans) <- 1
                                                                               ans <- t(sdev * ans) * sdev
                                                                           },
                                                                           factor={
                                                                               ans <- list(loadings=loadings, uniquenesses= 1 - comm,
                                                                                           sdev=sdev, call=match.call())
                                                                               class(ans) <- "statfacmodBurSt"
                                                                           })
                                                                        ans
    }
    
    # obtain covariance matrix
    # assign equal weight to all observations
    covmat <- factor.model.stat(data, weight=seq(1, 1, length=nrow(data)))

    # do factor analysis
    fa <- factanal(factors={NUMBER_OF_FACTORS}, covmat=covmat, rotation = "{ROTATION}")
    
    # since factanal doesn't do scores when given a covmat, do them manually
    # only for regression scores and orthogonal rotation?
    fa$scores <- scale(data, mean(data), sd(data)) %*% (solve(fa$correlation) %*% fa$loadings)
}

# write results to loadings.csv
try(write.csv(fa$loadings, quote=FALSE, file="{TEMP_DIR}/loadings.csv"), TRUE)

# write results to scores.csv
try(write.csv(fa$scores, quote=FALSE, file="{TEMP_DIR}/scores.csv"), TRUE)

